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I. Introduction 

In order to understand the global structure, 
dynamics, and physical and chemical 
processes occurring in the upper atmospheres, 
exospheres, and ionospheres of the Earth, the 
other planets, comets and planetary satellites 
and their interactions with their outer particles 
and fields environs, it is often necessary to 
address the fundamentally non-equilibrium 
aspects of the physical environment. These 
are regions where complex chemistry, 
energetics, and electromagnetic field 
influences are important. Traditional 
approaches are based largely on 
hydrodynamic or magnetohydrodynamic 
(MHD) formulations and are very important 
and highly useful. However, these methods 
often have limitations in rarefied physical 
regimes where the molecular collision rates 
and ion gyrofrequencies are small and where 
interactions with ionospheres and upper 
neutral atmospheres are important. 

At the University of Michigan we have an 
established base of experience and expertise 
in numerical simulations based on particle 
codes which address these physical regimes. 
The Principal Investigator, Dr. Michael 
Combi, has over 25 years of experience in the 
development of particle-kinetic and hybrid 
kinetic/hydrodynamics models and their 
direct use in data analysis as well as for 
theoretical studies. He has also directed and 
performed ground-based and space-based 
remote observational work and participated 
on spacecraft instrument teams. His research 
has involved studies of cometary atmospheres 
and ionospheres and their interaction with the 
solar wind, the neutral gas clouds escaping 
from Jupiter’s moon Io, the interaction of the 
atmospheres/ionospheres of Io and Europa 
with Jupiter’s corotating magnetosphere, as 
well as Earth’s ionosphere. 

Support from NASA's Applied 
Information Systems Research Program 
(AISRP) has resulted in the construction and 
first applications of three-dimensional, neutral 


and ion kinetic particle models, which have 
application to various space science problems. 

II. Progress and Current Status of the 
Kinetic Modeling Program 

As part of our AISRP project we have 
developed massively parallel code for 
performing kinetic calculations for neutral 
and ionized gases in space science 
applications. First we will describe some 
general work on the computational mesh. 
Then we will briefly describe an application 
to the plasma interaction of Jupiter's 
magnetized plasma torus with the satellite Io, 
and show some exciting new results. Thirdly, 
we will describe in more detail an application 
to the expanding neutral atmosphere of a 
comet. 

An Unstructured Computational Mesh 

The original plan at the beginning of our 
AISRP grant was to adopt the hierarchical 
Cartesian mesh algorithm used by our 
colleagues at the University of Michigan, lead 
by Dr. Tamas Gombosi. It became clear, in 
using and testing our own code, that the 
overhead in memory required for 
implementing their block-based parallel 
scheme, which works well for them in MHD, 
is far too inefficient in general for the more 
memory-intensive particle kinetic methods. 
At the same time, Dr. Gombosi's group itself 
has been in the process of interfacing their 
code with modules from other groups, by 
simply translating results between models 
where the spatial mesh of each model is 
optimized to the model in question. It was 
clear we should adopt the best computational 
mesh for our own problems. 

We spent considerable time developing an 
unstructured mesh, which is far more efficient 
for particle-based algorithms like ours. It is 
based on triangles for two-dimensional 
problems and tetrahedrons in three 
dimensions. This approach provides great 
flexibility in being able to conform to 
planetary-type problems where there is a 
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spherical inner body or an irregular comet as 
well as to very extended spatial regimes 
where a wide range of computational cell 
sizes is required. For example in Figure 1 we 
show a comparison of an old ( Combi 1996) 
geometrical-based structured radial-polar 2D- 
axi symmetric mesh used for an expanding 
cometary atmosphere with an unstructured 
triangular mesh. Because of the efficiency 
gained in parallelization and the large 
memory of modem computers, we now use 
far more computational cells than the PI used 
8 years ago. 






Figure 1. Comparison of a Structured and 
Unstructured Computational Mesh. 

The flexibility and power of this approach 
will be more apparent when we first apply it 
to a non-spherical central body. Comet 
Borrelly. See Figure 2. The triangles on the 
surface serve as the faces of the tetrahedrons 
at the boundary of the kinetic calculation. 
Science applications such as this, which are 
supported by NASA missions and science- 
topic specific Research and Analysis grants, 
illustrate the synergistic importance of the 
AISRP support for basic development of 
computer methodologies and algorithms. 


Figure 2. Image of Comet Borrelly from 
Deep Space 1. Below is the 3D mesh for an 
axisymmetric approximation to the shape of 
comet Borrelly. 

A Kinetic Simulation for the Interaction of 
Jupiter's Plasma Torus with Io 

For plasma particle simulations, we have 
been studying the interaction of Jupiter's 
plasma toms with the volcanic satellite Io. 
The kinetic particle method allows us to treat 
the important non-equilibrium aspects of the 
physics that are important in the near Io 
environment. From our first results it appears 
that the major effect is from the ring-beam 
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shape of the velocity distribution function for 
pickup ions near Io. Our best MHD 
simulation ( Combi et al. 2002) of the 
interaction compared with Galileo data is 
shown in Figure 3 and results from assuming 
that the ratio of the specific heats, p, is 2 
instead of 5/3. This at least accounts for fact 
that a gyration-dominated plasma should act 
as if there are only 2 kinetic degrees of 
freedom. 





Figure 3. MHD Model Results Compared 
with Galileo Measurements. From the top, 
down are plasma density, temperature and 
three components of magnetic field. 

In particle simulations the explicit 
modeling of the ion gyrations means that we 
can maintain ring-beam distribution of pickup 
ions along with their diamagnetic effect. 
Figure 4 shows the distribution of density, 
temperature and the N/S component of the 
magnetic field (Lipatov and Combi, 2004). 
We see that the diamagnetic effect of the 


pickup ions increases the temperature in the 
flanks of the wake (the two temperature 
peaks) and broadens the magnetic field 
perturbation, markedly improving even the 
best MHD result. Such anisotropic pressure 
effects are mediated by wave-particle 
scattering and collisional scattering with the 
neutral exosphere, which are modeled 
consistently and naturally in the kinetic 
simulation. 



- 6 - 1-20246 



X!cruM;i 


-*XJ | 



-* 4 -2 0 2 4 # 



4 -2 0 2 4 # 



MCHL. 

-6 4 -2 0 2 4 # 

XjlcnAi) 


Figure 4. Results of Particle Kinetic Model 
for the Galileo JO Flyby. From the top down 
are results of number density, temperature and 
the N/S component of the B field. The non- 
isotropic pressure effects of the gyration 
dominated pickup ions produce the 
temperature peaks and broad B field 
perturbation. 

A Kinetic Simulation for a Gaseous Coma 
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The cometary atmosphere is a unique 
phenomenon in the solar system. Owing to its 
negligible gravity, comets produce a highly 
variable and extensive atmosphere with a size 
much larger than the characteristic size of the 
cometary nucleus. As the comet approaches 
the sun, water vapor with other gases 
sublimates, generating an escaping 
atmosphere from the surface of nucleus, 
which consists of ice and refractory materials 
(rocky and organic solids and dust). 
Sublimating gas molecules undergo frequent 
collisions and photochemical processes in the 
near nucleus region. Except in the case of 
bright comets like Halley, Hyakutake, and 
Hale-Bopp on their close approaches to the 
sun, much of the time and for most of the 
coma, the gas production rate is small enough 
that the gas flow is in a high Knudsen number 
regime except for the immediate regions 
around active areas on the sunlit nucleus. 
[The Knudsen number is defined as the ratio 
between the collisional mean free path 
between molecules and the local length scale 
of interest in the problem.] This is especially 
true for the short period comets which are on 
the target lists for recent and upcoming 
NASA and NASA-participating ESA 
missions to comets. 

Two-dimensional axisymmetric models 
for six species (H 2 O, CO, OH, H 2 , O, H) have 
been performed using our new parallelized 
code. A hard sphere collision model was 
assumed for intermolecular collisions. For 
water-water collisions the viscosity equivalent 
cross section ( Crifo 1989) has been used. For 
other collision pairs, collision cross sections 
based on relevant atomic data ( Combi 1996) 
were used. 

We performed calculations based on 
conditions expected for the then target of the 
upcoming Rosetta mission to comet Wirtanen 
( Tenishev and Combi, 2002). The 
rescheduling to the sucessful launch this year 
to Comet 67P/Churyumov-Gerasimenko 
means it will have somewhat different 


specific characteristics, but the breakdown of 
simple pure fluid conditions will be otherwise 
similar to these results. 

The nucleus of Wirtanen was assumed to 
be spherical with a radius of 0.6 km and have 
axisymmetric gas emission. The gas flux from 
the dayside is controlled by absorption of 
solar radiation, and the flux varies as the 
square of the cosine of the subsolar angle. It 
was assumed that only two gaseous species 
sublimate from the surface: H 2 O and CO. 
Both species are in equilibrium with the 
surface, and the outflow mixture is partitioned 
with 80% for H 2 O and 20% for CO. As 
representative cases we performed 
calculations for the comet at its most active at 
perihelion of 1.05 AU (Q = 5 x lO^/s) as well 
as farther away at 1.5 AU (Q = 8 x 10 26 /s). 
For the nominal mission the comet spends 
most of the mission time at distances larger 
than 1.5 AU, and as we find, the collisional 
regime is already rarefied enough so that a 
particle approach is already needed. The 
larger and smaller production rate conditions 
are actually relevant for the new Rosetta 
target comet Churyumov-Gerasimenko at its 
perihelion distance of 1.5 AU and at 2.0 AU, 
respectively. 

In our paper (Tenishev and Combi, 2002) 
we describe in detail the comparison of the 
flow fields and collisional regimes for the two 
cases. The salient aspects can be told using 
only the Knudsen number plot for the 1.5 AU 
case as well as the radial velocity fields at 
both cases. These are shown in Figures 5 a, b, 
and c. What we find is that for the comet at 
1.5 AU, only the region near the nucleus and 
immediately above the active area is in a low 
Knudsen number (collisional fluid) regime. 
For the case at 1.05 AU (perihelion) the entire 
region near the nucleus is in a low Knudsen 
number regime and the flow conditions can be 
understood in terms of a hydrodynamic 
calculation. The effects of the two different 
collisional regimes are illustrated by the 
different radial velocity fields (Fig. 4 b and c). 
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The top plot shows the Knudsen number for the 
model at 1.5 AU where the shaded region is in 
the low Knudsen number (<0.1) state. The 
middle and bottom plots are contours of radial 
velocity for 1.5 and 1.05 AU, respectively. The 
differences are due to the rarefied density for the 
1.5 AU conditions. 

In the more rarefied case the adiabatic 
acceleration of the flow proceeds more 


slowly. The situation would be even worse 
for the comet at much larger heliocentric 
distances, where most of the Rosetta mission 
actually takes place. 


Dusty-Gas Hydrodynamics and Kinetics 

We already have incorporated multiple 
dust particle-size populations in the DSMC 
code. This has direct applications for the 
dusty-gas comae of comets and volcanic 
plumes on planetary satellites such as Jupiter's 
Io and Neptune's Triton. The dust particles 
can fragment and sublimate to serve as an 
extended source of gas. We have completed 
the testing of the dusty-gas kinetic calculation 
in the fluid limit, and it does reproduce 
published dusty-gas hydrodynamics 
calculations of Gombosi et al. (1986). The 
description of the Euler equations for a dust- 
gas cometary coma have been described in 
detail in many papers, including some by our 
own group at the University of Michigan 
( Combi et al. 1999, Korosmezey and 
Gombosi 1989). One can write down 
equation of conservation of mass, momentum 
and energy for the gas 
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respectively, where p is the gas mass density, 
u is the gas velocity, p is the gas pressure, v, 
and p, are the velocity and mass density for 


dust particles of radius a t . The term — is the 
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radius, a,. F is the gas-dust drag force, which 
is related to the forces on the individual 


particle size populations as , where 

>s 

the size dependent force is given by 
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the Boltzmann constant, T is the gas 
temperature and m is the gas mean molecular 
mass. The accommodation of gas via 
collisions with hotter dust yields the dust-gas 
heat exchange rate Q g d which is given as 

Q gd =^^ p u2(TriT i )St i 
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Here C p is the gas heat capacity at constant 
pressure and the rest of the coefficients can be 
defined under the assumption of diffusive 
reflection such that: 
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Finally, pd is the bulk mass density of dust 
particles of radius a,; T is the gas temperature, 
Ti is the dust temperature assumed in 
equilibrium with solar radiation. Otherwise 
the other intermediate quantities are s h the 
relative Mach number between gas and dust, 
Cd\ the dust-gas drag coefficient, T[ ec , the 
recovery temperature, R,\ the heat transfer 
function and St,', the Stanton number. Most 
of the dust-gas drag physics comes from the 
formulation by Finson and Probstein (1968) 
with later corrections discussed by Wallis 
(1982) and Kitamura (1986). 


We have already used the dusty-gas 
kinetic model for modeling the dusty gas flow 
relevant for the NASA Deep Space 1 mission 
flyby of comet 21P/Borrelly ( Tenishev and 
Combi, 2003). Results for the dust particle 
terminal velocities (i.e., the asymptotic 
velocities at large distances) as well as the 
acceleration profiles for different particle 



Figure 6. Results of Dusty-Gas Kinetic 
Calculations for Comet Borrelly during the 
Deep Space 1 Flyby. The upper plot shows 
the terminal velocity of dust particles as a 
function of dust particle sizes. The lower plot 
shows the velocities as a function of distance 
from the center of the nucleus. 
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III. Application of Developed Kinetic 
Modeling Tools for the NASA 
Clearly the model codes and algorithms 
developed here will serve as the core of our 
own theoretical modeling and data anaysis 
tools for now and future. The models are 
being used as part of the Pi’s ongoing NASA 
programs: (1) NAG5-12822 Studies 

ofTenuous Planetary Atmospheres, (2) 
NAG5- 13239 Distribution of Gas in the Inner 
Comae of Comets. These projects will 
continue to support science applications and 
publications of the code and algorithms 
developed. These were also used as the basis 
for 3 proposals written in 2004 and submitted 
to planetary R&A programs to the 2004 
ROSS NRA to study Jupiter's icy satellites 
and Mars' tenuous upper atmosphere, as well 
as for further proposed model and algorithm 
development under the AISRP. Model 
development and application is a 
sophisticated interdisciplinary and labor 
intensive effort, which is very difficult to 
support as minor parts of several focused 
space application oriented research grants. 
For the theoretical modeling community 
AISRP can provide something analogous to 
the Instrument Development programs (e.g., 
PIDDP), which are recognized by NASA as 
necessary to promote future basic 
development of instrumentation. 

Finally, one of the important goals of the 
AISRP is to provide new computational tools 
for use in future NASA programs. An 
interdisciplinary international modeling effort 
was begun this year, organized by the 
International Space Science Institute (ISSI) in 
Bern, Switzerland, in association with the 
Rosetta project (ESTEC) and the JPL Rosetta 
project (the support and organizing 
infrastructure for US investigators on the 
Rosetta project). As part of this effort, this PI 
is the lead investigator in the area of dusty- 
gas coma modeling (with T. Gombosi, U. 
Michigan and N. Thomas, U. Bern). The 
larger effort involves several other teams of 


investigators performing (1) MHD modeling 
of the comet/solar-wind interaction (T. 
Gombosi, K. Hansen, and M. Combi, U 
Michigan), (2) nucleus and surface-boundary 
modeling (B. Davidsson, ESTEC and H. 
Rickman), (3) electron energetics (T. Cravens 
U. Kansas), (4) plasma hybrid modeling (U. 
Motschmann and T. Bagdonat, MPI 
Braunschweig), (5) coma chemistry (K. 
Altwegg, U. Bern, and the entire U. Michigan 
group). Rosetta Project Scientist G. 
Schwehm and US Rosetta Project (JPL) 
Scientist C. Alexander are also participants as 
are a number of other Rosetta scientists. The 
dusty-gas DSMC model will be an integral 
part of this modeling effort for the next 12 
years or more. We have already begun a 
process which will lead to a set of coupled 3D 
models for the target comet in about a year, 
from the nucleus and out past the bow shock, 
including: nucleus porous upper surface 
layers, surface/coma Knudsen layer, dusty- 
gas coma, MHD and hybrid kinetic 
descriptions of the plasma interaction, coma 
chemistry, and electron energetics. 

Publications/presentations under AISRP 
Axisymmetric Dusty Gas-Kinetic Models for 
Comet 46P/Wirtanen. M.R. Combi, T.-M. 
Ho, and N. Thomas. American 
Astronomical Society, DPS meeting #33, 
#20.18, 2001. 

DSMC Simulation of the Cometary Coma. 
V.M. Tenishev and M.R. Combi. 
Proceedings of the International 
Symposium on Rarefied Gas Dynamics, 
Whistler, BC, July 20 - 25, 2002. 

DSMC Simulation of the Cometary Coma. 
M.R. Combi, and V.M. Tenishev 
Asteroids, Comets, Meteors, ESA 
Publications, Noordwijk, The 
Netherlands, 2002. 

3D Boltzmann Simulation of the Io's Plasma 
Environment with Adaptive Mesh and 
Particle Refinement. A. S. Lipatov and 
M.R. Combi. Eos Trans. AGU, 53(47), 
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Fall Meeting Suppl. Abstract #SM21A- 
0529, 2002. 

Development of a General Purpose 3D 
DSMC Flow Solver on Unstructured 
Meshes. V. Tenishev and M.R. Combi. 
AIAA Paper 2003-3776, 2003. 

Numerical Study of Dust and Gas 
Distributions in the Inner Coma of Comet 
Borrelly. V.S. Tenishev and M.R. Combi,. 
American Astronomical Society, DPS 
meeting #35, #38.01, 2003. 

3D Boltzmann Simulation of Io's Plasma 
Environment: Comparison with 

Observational Data. A. S. Lipatov and 
M.R. Combi. Eos Trans. AGU, 84(47), 
Fall Meeting Suppl. Abstract #P32A- 
1064, 2003. 

3D Hybrid Simulation of the Cometary 
Plasma Environment with Adaptive Mesh 
and Particle Refinement. A. S. Lipatov 
and M.R. Combi. Eos Trans. AGU, 
84(47), Fall Meeting Suppl. Abstract 
#SM31C-1121, 2003. 

DSMC Simulation of the Cometary Coma. 
V.M. Tenishev and M.R. Combi. Rarefied 
Gas Dynamics: 23rd International 
Symposium. AIP Conference Proceedings, 
663, 696-703, 2003. 

3D Boltzmann Simulation of the Io's Plasma 
Environment with Adaptive Mesh and 
Particle Refinement. A. S. Lipatov and 
M.R. Combi. Icarus (submitted), 2004. 

Effects of kinetic processes in shaping Io's 
Global Environment: 3D hybrid and fluid- 
ion-particle-ion models of the Galileo 127 
Flyby. A. S. Lipatov and M.R. Combi. 
Icarus (in preparation). 
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Development of a General Purpose DSMC Flow Solver on Unstructured 

Meshes (AIAA Paper 2003-3776, 2003) 

V. Tenishev*, M. Combi 1 

Space Physics Research Laboratory, Department of Atmospheric, Oceanic and Space Sciences, University 

of Michigan, Ann Arbor , MI 


Abstract 

A general purpose kinetic parallel flow 
solver was developed based on DSMC 
methodology. Basic characteristics of the solver 
and some examples of its application are outlined 
in the paper. 


Introduction 


One of the most important aspects of the 
rarefied gas simulation is a high value of the 

Knudsen number Kn = Aj L in the flow and, 

consequently, the requirement of solution of the 
Boltzmann equation 


I • • 




The situation is usually more complicated 
because most of the practical applications 
involve consideration of diatomic gases, which 
can be excited vibrationally and rotationally and 
for which ionization, dissociation and chemical 
transformations are possible. A need to model 
non-equilibrium hypersonic flows motivated the 
development 1 of numerical techniques to treat 
such problems. Nowadays the DSMC method 2 is 
de facto the standard method for rarefied gas 
dynamics, where the state of the rarefied gas 
flow, which is presented by an assembly of 
molecules, is determined by collisional dynamics 
of a finite number of model particles and, hence, 
hold potential for providing information on flows 
where the collisional rate is not sufficient to 


maintain equilibrium energy distribution. 


* Research Assistant; vtenishe@umich.edu 
f Senior Research Scientist 


This paper presents the current status of 
development of a general purpose DSMC flow 
solver and some results of its application. 

Code 


Current problems of gas dynamics 
cannot be solved by methods other than 
numerical ones. A wide range of practically 
important problems requires developing a highly 
general numerical tool. C++ language was 
chosen in the present work. 

Requirements 


The following main principles guided the 
development of the code: 

1 . Portability: the code was tested on a set 
of different platforms and compilers. 

2. Easy modification of the core of the 
code: change of the basic physico- 
chemical models installed into the code 
should not require any changes in any 
other part of the code. 

3. Careful usage of resources: special 
attention has been paid to optimize 
allocation of memory and CPU 
resources. 

4. Parallel implementation: MPI library 
has been used in parallelization of the 
system. Dynamic and static load 
balancing capabilities were 
implemented. 

5. Easy adaptation of the code to particular 
physical problem (without modifacation 
of the core of the code). 

6. Friendly user interface: easily 

understandable input and output files. 

7. Pre/post-processing: including pre/post- 
processor into the core of the system. 
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allow running and post-processing on 
different platforms. 

8. Optimization of the numerical 
procedure for the current problem: 
change main parameters of simulation 
on order to improve overall 
performance. 

Basic Specification of the System 

The code potentially can be applied to a 
wide range of problems. The following is the list 
of its main features. 

1. Neutral gas flow analysis in 1, 2, 3D 
cases. Spherical- and axis-symmetry 
can be taken into account. 

2. The flow field is reconstructed on an 
unstructured mesh. 

3. Multiple species gas flows can be 
modeled. Chemical reactions and 
relaxation of internal degrees of 
freedom can be considered. 

4. A set of molecular and collisional 
models is available. 

5. Mean free path and model particle 
weight adaptation has been 
implemented 

Code Adaptation for Particular 
Problem 

The core contains a general set of 
models (molecular models, collisional, chemical 
models). Some cases require introduction of 
additional particles or/and physicochemical 
models to solve particular problem. Such models 
may contain a set of approximations or/and 
experimental data specifically suited for the 
problem and, so, they may not be useful for any 
other cases. It is good to have the possibility to 
incorporate models into numerical procedure 
without any modification of the core of system. 
There are several difficulties that have to be 
solved in this way. Based on capabilities that are 
provided by C++, an approach of treating 
problem-dependent models was developed and 
incorporated into the system. 

Optimization techniques 

The accuracy of the numerical results, which 
are obtained by the DSMC method, is very 
dependent on the size of the sample and, finely, 
on the number of model particles per cell. So, in 
order to make the most effective use of computer 


resources, it is desirable to get the distribution of 
model particles over the domain close to 
uniform. This problem is especially important 
for the case of a 2D axi symmetrical simulation. 
There are several methods 2 to control the 
distribution of particles: variation of time steps, 
grid manipulation and direct variation of a 
particle’s weight. 

Several techniques could be used to decrease 
total computational cost of steady state 
simulation. Due to methodology of the DSMC 
method, a steady state solution can be obtained 
as the limit of unsteady flow simulation, where 
convergence takes a considerable part of total 
computation time. To get faster convergence, it 
is possible, starting with low accuracy, to 
improve it gradually as the current solution 
converges to the finial result. One of the possible 
implementations 3 of the strategy is to gradually 
increase total number of model particles as the 
flow pattern converges to the steady state. 
Another scheme was developed and 
implemented in the current work. The simulation 
is conducted on a set of meshes with different 
levels of refinement, starting with a coarse grid 
at the very beginning and maintaining the 
number of model particles per cell in desired 
limits. As the number of collisional cells and 
particles employed increase, the resolution of the 
computational procedure improves both in terms 
of flow macroparameters distribution and 
captured physical processes. 

The use of parallel computers can 
significantly increase the range of application of 
the DSMC method. Effective use 4 of such 
computational systems can be achieved only 
with careful load balancing. 

Static and dynamic load balancing 
capabilities are incorporates into the current 
implementation. Different criteria are available 
for the balancing. Static load balancing may be 
based either entirely on parameters of initial flow 
condition or on some knowledge about the 
structure of converged solution. In the first case, 
the load is balanced by equipartitioning of total 
number of particles, total volume or total cell 
number among all available processors. In 
addition to them, a total execution time can be 
used as a criterion in dynamic load balancing. 


2 


Applications of the code 

Cometary Comae Simulation: Generic 
Jupiter Family Comet 

The cometary atmosphere is a unique 
phenomenon in the solar system. It represents 
highly extended dusty gas cloud with a large 
range of change of its macroparameters from 
fluid to collisionless. So, kinetic methods are the 
most suitable tools for simulation processes in 
cometary comae. 

PHYSICAL MODEL 

When the nucleus approaches the sun, it 
absorbs solar radiation, resulting in an increase 
of gas sublimation from the surface. Because the 
gravitational force is negligible, the gases 
leaving the surface form a cometary exosphere. 
The radiation reaching the surface and supplying 
energy for sublimation is partially adsorbed 
dusty atmosphere. Any changes in the gas and 
dust production result in changes of optical 
characteristics of the cloud. The main volatile 
component of the cometary coma is water. So, 
water sublimation controls sublimation of others 
components when a comet is within about 4 
astronomical units (AU) of the sun. Freshly 
evaporated molecules are photodissociated and 
photoionized, and therefore most of the chemical 
kinetics of cometary atmospheres involves the 
resulting highly reacting ions and radicals 5 . Most 
of the absorption of solar radiation takes place at 
the surface or in the dense region near the 
surface. Models of temperature and gas 
production distributions on the surface of a 
spherical nucleus can be found 6,7 . The simplest 
model of gas production assumes an optically 
thin coma, which is an excellent approximation 
for most comets, most of the time. 

The conditions for a comet coma 
usually consider the gas at the surfaces to be 
stationary at a given temperature. This 
assumption neglects the real physics of the gas 
surface interface. But to get correct structure of 
the flow field in the near nucleus region, it is 
more important that the flux integrated around 
the surface equals the desired production rate and 
the gas temperature near the surface is close to 
vaporization temperature of water. Due to 
absence of gravity, the downward flux is 
negligible. This result in upward-half 
Maxwellian distribution of water vapor in a gas 
layer where thermodynamic equilibrium is not 


held. So, in a cometary atmosphere, unlike that 
of planets atmospheres, the mean flow velocity 
can be comparable with to the thermal speed 5 , 
therefore methods adopted for studying planetary 
atmospheres are not longer valid. The flow then 
quickly thermalizes and passes smoothly through 
the sonic transition. 

Two-dimensional axi symmetrical 
models for six species (H 2 0, CO, OH, H 2 , O, H) 
were considered in this work. A hard sphere 
collision model was assumed for intermolecular 
collisions. For water-water collisions the 
viscosity equivalent cross section 8 has been used. 
For other collision pairs, collision cross sections 
based on relevant atomic data were used 9 . 

The cometary nucleus was assumed to 
be spherical with a given radius and 
axisymmetrical gas emission. The gas flux from 
the dayside is controlled by absorption of solar 
radiation, and the flux is varying as the square of 
the cosine of the subsolar angle. It was assumed 
that only two gaseous species sublimate from the 
surface: H 2 0 and CO. Both species are in 
equilibrium with the surface and in the outflow 
mixture is partitioned on 80% for H 2 0 and 20% 
for CO. 

Photochemistry is included into the 
model. The main reaction branches are presented 
in Table 1 


Reaction 

Branching 

ratio 

h 2 o-*h+oh 

0.88 

h 2 o^h 2 +o 

0.22 

OH — » O+H 

1.00 

H 2 H+H 

1.00 


Table 1 Photochemical branching 

STATEMENT OF THE PROBLEM 

A comet with a production rate of 

3xl0 28 molecules/sec was chosen, as typical of 
Jupiter family comets which will be spacecraft 
mission targets over the next decade or two. The 
radius of the nucleus was taken to be 10 km and 
the gas production rate was distributed as the 
square of the cosine of the subsolar angle with 
no gas emission on the night side. The gas was 
initiated at the surface at rest and with a 
temperature of 180K. The results of the 
calculation are shown in Fig. 1,2. 

For photochemical reactions, the rates 
are typically much smaller than collision rates 
near the nucleus. Taking the characteristic time 
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constant for particular photochemical reaction T. 
the probability of photodecay for a single 
molecule at a given time step A t is 

Pi =l-exp(-A tfc) 

The regime of expansion of the 
cometary coma changes from collision 
dominated near the surface of the nucleus though 
transitional one to collisionless in the far field. 
Consequently, all of the regimes have to be 
described by numerical approach chosen for 
simulation, which makes questionable 
applicability of fluid based approximations. 



Fig. 1 The two-dimensional field of Knudsen numbers for the 
case of a generic Jupiter family comet. 



Fig. 2 The two-dimensional field of decimal logarithm of 

water density (in cm ) for the case of a generic Jupiter 
family comet. 


Gas/Dust Interaction in cometary 
comae 


Our present understanding of cometary 
nucleus is based on idea of “dirty snowball” 
according to which the nucleus consists of frozen 
volatiles and nonvolatile dust. The sublimated 


dust molecules drag away dust particles from the 
surface of cometary nuclei, so, the gas/dust 
interaction may a have significant effect on 
macroparameter distributions in comeraty 
comae. In this work, our particular interest was 
concentrated on modeling of inner comae 
processes, the region where dust and gas are 
produced and accelerated to their terminal 
velocities. 

The following mechanism of dust 
particle acceleration is adopted here. It is 
assumed that an influence of gravitational force 
is negligible and, so, dust particles acceleration is 
only due to collisions with the gas molecules 

Different kinds of gas molecules/dust 
surface interaction mechanisms can be easily 
introduced into the model, the simplest version 
of which assumes elastic reflection on the 
surface of the dust particle. In the current work, 
its simplified approximation has been used. 

Taking a as the dust particle radius, V dust and 

V mo/ as dust and gas particle’s velocities, the 
change in momentum of the dust particle during 
time interval tSX can be expressed as 
= 

= A |V to; ~V dust )^ 

mol ''cell 

Here, F mol is the model particle weight of gas 
molecule and V cell the volume of the 

computational cell. The model has been 
validated via comparison with test solution 
obtained by a fluid/dust model 5 . For the 
problems of cometary comae simulation, due to a 
wide range of collisional regimes in the flow, 
application of the fluid approximation is 
restricted. 

The problem was considered in ID 
spherically symmetrical geometry. The nucleus 
radius R = 5. 3km and gas production rate of 

3.5 XlO 28 molecules/s corresponds to 
parameters of comet Borrelly at heliocentric 
distance of 1.36 AU. It was assumed that only 

H 2 0 molecules are sublimates from the surface. 

The dust/gas mass ratio is 0.4. 

Some obtained results are shown in Fig. 

3, 4. 
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of its radius. Here U w is the terminal velocity of the gas 
phase. 



Fig. 4 Radial mean velocity (m/s) profiles for dust particles 
of a given size and the gas phase in the cometary coma as a 
function of a distance from the nucleus. 


3D Two Body Interaction Problems 

The ultimate goal of the work is 
developing of a general purpose DSMC flow 
solver that can be used to solve a wide range of 
problems in complex geometries. A two-body 
interaction problem has been chosen to 
demonstrate the capability of the code to treat 
complex 3D flow problems. The geometry of the 
problem is presented on the Fig. 5. 



Fig. 5 Geometry sketch for the 3D two body interaction 
problem 


Two bodies of revolution are 
cone/cylinder combination. The cone semi-angle 

is P — 20° , the cylinder length is 

L = 250mm , its diameter is d — 50 mm . The 
distance between the body axes is 
A z = 70mm . The bodies of revolution were 
aligned parallel, without relative displacement at 
a zero incidence angle. 

The flow of N 2 has been simulated on 
a mesh with 5.7 X 10 6 computational cells using 
— 70xl0 6 model particles. The upstream 
conditions are, U„ = 1 242 .6m / s , 

T„=51K, n„= 8.95X10 21 , which 

corresponds to Mach number of M=9.91. Some 
of the results are presents of Fig. 6, 7. 



Fig.6 Normalized number density profile at a distance of 
0.61 L from the leading edge of the body 
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Fig. 7 Normalized temperature profile at a distance of 0.67L 
from the leading edge of the body 


Conclusion 

The current state of developing a 
general-purpose kinetic flow solver was 
presented. The tool has been validated via 
comparison with solutions of a set of test 
problems. 

Main features of the code were outlined. 
It was stressed that portability, modularity, easy 
adaptation to particular problem, parallel 
implementation and friendly user interface 
makes the developed code a useful tool for 
analysis of complex multidimensional gas flows 
of chemically reactive mixtures. Discretization 
of the flow field on unstructured meshes allows 
the code to be used for solution of problems of 
internal and external gas flow in the case of 
complex geometries. Adaptation capabilities can 
be utilized to optimize overall performance. 

In the paper, some examples of 
application of the code have been presented. A 
problem of two-phase dust/gas flow has been 
considered. 
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Abstract. The global dynamics of the ionized and neutral components in the 
environment of Io plays an important role in the interaction of Jupiter’s corotating 
magnetospheric plasma with Io. The stationary simulation of this problem was done in 
the MHD and the electrodynamics approaches. One of the main significant results 
from the simplified two-fluid model simulations (Saur et al., 2002) was a 
production of the structure of the double-peak in the magnetic field signature 
of the 10 flyby that could not be explained by standard MHD models. In this 
paper, we develop a method of kinetic ion simulation. This method employs 
the fluid description for electrons and neutrals whereas for ions multilevel, drift-kinetic 
and particle, approaches are used. We also take into account charge-exchange 
and photoionization processes. Our model provides much more accurate 
description for ion dynamics and allows us to take into account the realistic 
anisotropic ion distribution that cannot be done in fluid simulations. The first 
results of such simulation of the dynamics of ions in the Io’s environment are discussed 
in this paper. 

KEYWORDS: Io; Jupiter; Magnetosphere; Satellites; Atmospheres. 
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1. Introduction 

The global interaction of the Jovian plasma torus with the Io is a fundamental 
problem in magnetospheric physics. It requires the solution of a highly nonlinear 
coupled set of integro-MHD/kinetic-Boltzmann equations which describe the dynamics 
of Jupiter’s corotating magnetospheric plasma, pickup ions and the ionosphere together 
with the atoms from Io’s environment. To leading order, the plasma and neutral atoms 
are coupled by resonant charge exchange, although other coupling processes are present. 
The characteristic scale of the ionized components is determined usually by the typical 
ion gyroradius, which for Io is much less than characteristic global magnetospheric scales 
of interest. By contrast, the mean free path of neutral particles and the typical ion 
gyroradius are comparable to characteristic exospheric scales such as the thickness of 
the exosphere and atmosphere, the distances separating the possible magnetic barriers 
and the surface of Io, etc. Consequently, the Knudsen number Kn = \/L exo (A, the 
mean free path of neutral particles and L exo , a characteristic exospheric scale), which is 
a measure of the distribution relaxation distance, satisfies Kn as 1. Thus it is difficult 
to assume that the distribution of the neutral atoms and ions can relax to a Maxwellian 
distribution, and one needs ideally to solve a Boltzmann equation for the neutral 
and ionized components in which charge exchange and photoionization processes are 
included. 

Several approaches for including the neutral component and pickup ions self- 
consistently in models describing the interaction of plasma torus with Io have been 
formulated. Early theoretical work was often done either in the context of a “thin” 
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atmosphere (e.g., see Cloutier et al., 1978) indicative of the surface temperature ( 130 K), 
or “thick” extended neutral atmosphere (e.g., see Goertz, 1980) more indicative of 
volcanic temperatures ( 1000K). Subsequent evidence (see the review by Lellouch, 1996) 
seems to indicate a mixed picture of the global atmosphere, which has a large extended 
corona like a thick atmosphere, but appears to be dominated by local major injection of 
hot (high speed) gas/dust plumes to high altitudes but only near active volcanic vents. 
Therefore, although the atmosphere is probably only locally thick, it still has a large 
extended neutral corona which might provide a sufficient source of impact ionization 
and photoionization to explain the plasma torus. 

Southwood et al. (1980) examined data from several Voyager instruments and 
examined the possible role of an intrinsic magnetic field for Io as a way to retain a 
robust enough ionosphere, which could provide enough conductivity for completing the 
Io-Jupiter current circuit. Neubauer (1980) presented an analytical model of 
the Alfven standing wave current system which connects current through 
the ionosphere of Io. Southwood and Dunlop (1984) and Ip (1990) suggested models 
stating that mass loading effects should result in the formation of a tail-like structure in 
the wake behind Io; as a consequence of the enhanced plasma density, the magnetic field 
perturbations are continued into the wake. Thus the field-aligned current is not only 
generated by Io itself, but also in the wake. Several years after Voyager, 3D numerical 
studies of the plasma flow past Io were performed using electrodynamic (Wolf-Gladrow 
et al., 1987), magnetohydrodynamic (Linker et al. 1989; Linker et al. 1991) and resistive 
magnetohydrodynamic (Kopp, 1996; Kopp, Birk and Otto, 1998) approaches. 
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There have been recent efforts to improve and extend the pre-Galileo simulations 
both in terms of the MHD (Combi et al., 1998; Linker et al., 1998; Kabin et ah, 2001) 
and the electrodynamic (Saur et al., 1999) approaches. These two approaches are 
distinguished by the physical assumptions which they each do and do not (or in some 
cases, can and cannot) include. However, MHD cannot, at least yet, include the effects 
of realistic conductivities (Hall and Pederson) or charge separation effects which are 
likely to be important very close to Io where the neutral densities are large and electric 
potential can introduce non-symmetric flow around body. They either include constant 
artificial conductivity (Linker et., 1998) or assume perfect conductivity (Combi et al., 
1998), however comparisons of the sets of published results do not indicate that this 
choice has any important consequence. 

The non-magnetic models produce magnetic field perturbations that are similar to 
the Galileo measurements, but none are quite as deep or as broad, and none have the 
reversal of the perturbation (the double-peaked structure) in the center of the wake. The 
magnetized models of Linker et al. (1998) produced a broad and deep perturbation, but 
not the self-reversal at the center of the wake (the double-peak or bite-out). Thus, the 
observational data, e.g., density and magnetic field profile along the Galileo trajectory 
cannot be fully explained by MHD models. But the most recent evidence from all the 
Galileo flyby data is that Io does not possess any substantial internal magnetic field of 
its own (Kivelson et al., 2001). 

In the electrodynamical model of Saur et al. (2002), they added an extra ionization 
source, based on the high-energy bidirectional electrons observed by Williams et al. 
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(1996,1999) and Frank and Paterson (1999). These high-energy electrons were included 
as an energy source in addition to electron impact ionization by the thermal electrons 
and photoionization to create a dense plasma wake as is observed. Simulations without 
these extra bidirectional electrons show a nearly empty wake Saur et al. (1999). They 
also concluded that the electrodynamic part of Io’s interaction is best described as an 
ionosphere-like interaction rather than a comet-like interaction (Saur et al., 2003). Saur 
et al. (2002) demonstrated first that the diamagnetic and inertia currents 
are responsible for formation of the double-peak magnetic profile along the 
10 trajectory. They received the magnetic field profile that has an oscillating 
structure and the maximum in the magnetic field profile is to some extend 
was a narrower than in observation for standard model of the atmosphere 
(Fig. 11, Saur et al., 2002). In a case of longitudinally symmetric atmosphere 
they received a strong double-pick structure but the spatial scale of this 
structure was a twice smaller than in observation (Fig. 14, Saur et al., 2002) 
and the space scale of these peaks are much smaller than in observation. 

In this paper, we describe a new approach to solving the time-dependent Boltzmann 
equation together with a hybrid plasma (ion kinetic) model in three spatial dimensions 
using a prescribed atmosphere model for Io. A Boltzmann simulation is applied to model 
charge exchange between (incoming and pickup) ions and immobile atmosphere. Several 
simulations are run, and the results described. We show for the first time the predicted 
distribution for ions and the electromagnetic field throughout Io’s environment. The 
results of such kinetic simulations are compared with those obtained from a related 
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MHD model and observational data. We also found the best parameters for our model 
that describe satisfactory the main behavior of the observational data at least on the 
qualitative level. 

2. Formulation of the Problem and Mathematical Model 

2.1. Simulation Model 

To study the interaction of the plasma torus with the ionized and neutral 
components of Io’s environment we use quasineutral hybrid models for ions and 
electrons. The interaction of the ions with atmosphere is dominated by charge exchange. 
The atmosphere is considered to be an immobile component in this paper. The general 
scheme of the global interaction of the plasma torus with Io and the Galileo 10 trajectory 
is given in Fig. 1. The Galileo 10 flyby occurred nearly in the equator plane of Io and 
perpendicular to the direction of the plasma wake defined by the corotating plasma flow 
past Io. The spacecraft trajectory passed approximately 900 km down-stream of Io (in 
the sense of the plasma torus flow). In our coordinate system the Z axis is parallel to 
Uq (corotational plasma velocity), Y is aligned with the spin axis, and X = Y x Z. The 
relation between our coordinate system and IphiO (X*,Y*,Z*) system (see Fig. 3 from 
Kivelson et al., 2001) is the following: X is parallel to Y*, Y is aligned with Z*, and Z 
is parallel to X* . 

In the hybrid simulations described here, the dynamics of upstream ions and 
implanted ions is described in a kinetic approach, while the dynamics of the electrons is 
described in a hydrodynamical approximation. 
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The single particle ion distribution function / s (t,x,v) has to fulfill the Vlasov 
equation 

d , d , F5, 

where F symbolizes all forces acting on the ions. 

The single ion particle motion is described by the equations 


d r 


S,l 


= v S) b 


dv 3i , _ Z a e / v 3 j x B\ 
dt M s V + c ) ’ 


E* = E - 


( 2 ) 


dt dt 

Here and M„ denote the charge state and the mass of the ions, the subscript 
s denotes the ion species (s = 1 for incoming ions and s = 2 for pickup ions) and the 
index l is the particle index. a e ff is an effective conductivity that may include Coulomb 
collisions and collisions due to particle-wave interaction, in the sense of a standard 
resistive Ohm’s law. Note that the resistivity used in Eq. (2) must depend on individual 
velocities of ions and electrons. However, we use the effective resistivity and an effective 
electric field E* in Eq. (2) to satisfy momentum conservation for the plasma system 
(Zueva et al., 1975). 

In the nonradiative limit Ampere’s law is given by 

47T 


V x B; 


( 3 ) 


and the induction equation (Faraday’s law) by 

15B 


c dt 


+ V x E = 0. 


( 4 ) 


The total current is given by 


J = J e + Ji; Ji = 53 Z s en s \J s , 

S = 1 


( 5 ) 
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here U s is the bulk velocity of ions of the type s. 

We further assume quasi-neutrality 

2 

n e = Yl Z s n s- ( 6 ) 

5=1 

For massless electrons the equation of motion of the electron fluid takes the form of 
standard generalized Ohm’s law (e.g. Braginskii, 1965): 

E = — (J e X B) - — Vp e + —Re, (7) 

en e c en e en e 

where p e = nm e (v ^)/ 3 = n e k B T e , and v' e are the scalar electron pressure and the thermal 
velocity of electrons, and the electron current is estimated from Eq. ( 5). 

The term including R e , symbolizes the mean momentum change of the electrons 
due to their collisions with the ions, which can be parametrized in terms of the collision 
frequency u es between electrons and ions of the species s as 

Re — -n e m e u es (U e - U s ) (8) 

Since we suppose that electron heating due to collisions with ions is very small, the 
electron fluid is considered adiabatic: 

Pe oc n 5 J 3 . (9) 

The ion kinetic approach allows us to take into account the effects of anisotropy of ion 
pressure, the correct mass loading processes, a penetration of ions across the ionosphere, 
and the asymmetry of plasma flow around the Io. Remember that the fluid models 
which account only for the scalar ion pressure may result in an extra- effusion of the 


pickup ions along the Alfven wing. 
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Charge exchange may be included in a hybrid model by a simple procedure (Lipatov 
et al., 1998; Lipatov, 2002). The total loss rate in s _1 for the ionized component has the 
form 

&x(x, V, t) = J / 0 (x, V tt , t)V re l, a a e x{Vrel,a)d 3 V a , (10) 

where the velocity of ions relative to an atom with velocity v a is V re ^ a — |v — v a |. 

If we suppose that the neutral component has a Maxwellian distribution, then {3 ex 
may be approximated as (Ripkin and Fahr, 1983) 

0ex(x, v, t) ~ n 0 (x, t)Vr e l, a <J ex (Vrel,a)- (11) 

In this paper we assume that the neutral component of Io’s atmosphere has zero bulk 
and thermal velocities. If one needs to include the nonzero temperature and bulk 
velocity, the effective average velocity of ions relative to an atom with velocity may be 
found from the general equation of Ripkin and Fahr (1983). 

Let the time interval t* with respect to the charge exchange be a random variable 
with a distribution function 

wi(t*) = exp[— f Pexjdt], (12) 

Jto 

where / is the particle index. 

A survival probability against the charge exchange event w ex may be written as 
follows: 

Wex,l(t) = ex P[~ / Pex-Jtdt]. (13) 

Jto 

The integration is over the trajectory of the particle with index l. At the time of 
creation (either at the boundary of the calculation box, or at the moment of charge 
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exchange), an ionized particle has initial coordinates x;(t 0 ) = x ij0 , V((t 0 ) = v; ]0 , a weight 
a t (to ) = aifl, and a survival probability w ex j(t Q ) = 1. For each new ion with index l, we 
have to determine the critical probability w* ex l when charge exchange will occur, and 
this is done using the relation 

<,/ = 6 ( 14 ) 


where £ is random number on the interval from 0 to 1. During the calculation we have 
to identify those particles for which the probability of survival satisfies the condition 


^ex,i < W* ex l . (15) 

If the particle satisfies the condition (15), then we have to exchange the velocity of this 
ion with the velocity of an atom from the atmosphere of Io. In the present simulations, 
we do not take the cross section a into account in (10) because of the weak dependence 

on \v — v p \ (as was done in the work of Malama (1991)). If charge exchange occurs, then 

a new neutral particle begins its motion with v* = 0 and w eXf t = 1. 

Io’s environment model includes 2 sources of pickup ions: the extended coma halo 
distribution with pickup ion production law 

G OC ^i^atmosbke^r (lb) 


for 1 x ri 0 < r < 7 x r Io and the exosphere distribution close to Io with ion production 
law 

G oc t'in atmos W int exp[— (r - r Xo ) / H^ mos }. (17) 


Here W ext and W int denote the fraction of the total ion production rate in the halo 
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and exosphere (W tI , t — (100 — 95)% of the total), respectively. n atmos denotes the 
maximum atmosphere density, is effective ionization rate. 

The inner region of the Io’s ionosphere is described by immobile ions with the 
following density distributions: 

0.5n iono exp(— (r - r lo )/H efS ), for r > r Io , 

rciono(l - 0.5exp(-|r - r Io |/// e //)), for r < r io , 

We assume that the incoming flow has a small finite resistivity to suppress the “shot” 
noise fluctuations. We also take into account the effect of finite conductivity of Io’s 
body so that 

<7e// = CTup, for r > rio, 

°eff = CTio, for r < rio, 

Our code solves equations (1) - (8), (5) - (9), (10) - (15) and (16) - (17). 

Initially the computational domain contains only supersonic plasma torus flow 
with a homogeneous spatial distribution and a Maxwellian velocity distribution; the 
pickup ions have a weak density and spherical spatial distribution. The magnetic 
field and electric fields are B = Bo and E = Uo x Bo. Inside the Io electromagnetic 
fields are E =0 and B = Bo- In the cases examined in Figs. 9-16 we choose 
B 0 /Bo = (0,1,0); Hjo/Eq = (1,0,0) and in the cases examined in Figs. 2-8, 17-19 we 
choose B 0 /B 0 = (0.0394,-0.9854,0.1657) and E 0 /E 0 = (-0.98,-0.0394,0.0). 

At t > 0 we begin to inject the pickup ions with a distribution according to Eqs. 
(16-17). Far upstream (z = —DZ/ 2), the ion flux is assumed to have a Maxwellian 
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: distribution, 

i / = noo(7ru t 2 h )" 3/2 exp -- V ^ , (18) 

lV th . 

where v t h and U are the thermal and the bulk velocities of the plasma torus flow, 

respectively. 

i 

Far downstream, we adopt a free escape condition for particles and Sommerfeld’s 
radiation condition for the magnetic field. On the side boundaries (y = ±DY/2 and 

| 

x = ±DX/ 2), periodic boundary conditions are imposed for incoming flow particles and 
the electromagnetic field. In some cases we also tested the use of the upstream boundary 
condition for electromagnetic field on the side boundaries. In these situations we also 
employ a buffer zone with thickness about of 10 x Ax where a smoothing procedure 
provides a transition for electromagnetic fields from the perturbed value to the upstream 
value on the side boundaries (see e.g. Umeda et al., 2001), effectively allowing Io- 
generated Alfven disturbances to propagate away. The pickup ions come out from the 
computational domain when they intersect the surfaces x = 5 x Ax, x — DX — 5 x Ax, 
y = 5 x A y, y = DY — 5 x Ay, z = 5 x A z, z — DZ — 5 x A z. Thus the flux of 
the pickup ions is absent on the side boundaries. At Io’s surface the particles may be 
reflected or absorbed. Note that the position of the Io is x = 0, y = 0, z = 0. 

The three-dimensional computational domain has the dimensions DX = 20 L, 

DY = 20L, and DZ — 10L, or DX — 14 L, DY = 14 L, and DZ = 12 L, or DX = 10 L, 
DY = 10 L, and DZ = 8 L, where L equals the radius of Io, r\ Q = 1800 km. We use the 
meshes of 201 x 201 x 101, or 161 x 161 x 121, or 141 x 141 x 121 grid points, and 2 x 10 8 
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and 5 x 10 7 particles for ions and pickup ions, respectively, for a homogeneous mesh 
computation. The time step At satisfies the condition v max At < min(Ax, Ay, Az)/8. 

The relationship between dimensional (U, E , B , p e , n, T) and dimensionless ({/', 
E ! , B f , Pg, n 7 , T f ) parameters may be expressed via dimensional upstream values as 
follows: 

U = UVo, E = E'B 0 Uo/c, B = B'Bo, Pe = Pe'Pe o, 

n = n'n 0 , T = T'MJJl , (19) 

whereas the dimensional time and distance may be expressed via the bulk velocity f7o 
and characteristic scale L: 

t = t'L/Uo , ^ (20) 


2.2. Numerical Method 

We employ a standard particle-in-cell (PIC) method in the case of a homogeneous 
grid. The time integration of the particle motion equations uses a leapfrog scheme. The 
time integration of the electromagnetic field equations uses an implicit finite difference 
scheme (see, e.g., Lipatov (2002)). We use different time steps for particle and field 
pushing (subcycling). This code was optimized for parallel computation using MPI and 
OMP. 

Since the gyroradii must be resolved, a grid point spacing of less than 1 gyroradius 
is required in order to avoid numerical dispersion and dissipation. On the other hand, 
good statistics are required, therefore a sufficiently large number of particles per cell 
have to be used (i.e. low “shot” noise). We use a homogeneous mesh for the simulations 
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presented here. 

3. Results of the Simulation 

To study the interaction of the plasma torus with the ionosphere of Io the following 
set of Jovian plasma torus and ionosphere parameters were adopted in accordance 
with observational data (Frank et al., 1996; Kivelson et al., 1996): corotation velocity, 
Uq = 56.8 km/s; plasma density, no = 3500 cm -3 ; plasma temperature, 92 eV; mean ion 
mass, Mi = 22 M p ; ratio of specific heats, 7 = 5/3; ion and electron betas, $ = 0.039; 

= 0.0022; magnetic field, B = 1800 nT; Alfven and sonic Mach numbers Ma = 0.4 
and Ms = 2.2. 

For pickup ions we use the following parameters: total ion production rate, 

Q ion = (0.7 — 2.) x 10 28 s -1 ; mean ion mass, Mpi = 22 M p ; electron exosphere and 
ionosphere betas, /3 e ,ex 0 = 0 — 0.0022, /3 e ,iono = 0 - 0.0005; effective cross section for 
charge exchange, a exc h = 1-5 x 10 -15 cm 2 ; effective ionization rate, 17 = 10 -5 s -1 ; 
atmosphere scale height, H^ mos = 0.06 — 0.09 x ri 0 = 108 km— 162 km; maximum value 
of the density of the atmosphere, n atmos = (0.5 — 10) x 10 8 cm -3 ; effective ionosphere 
scale height, H e ff = (0.0001 — 0.00025) x r\ Q =(0.18-0.45) km. Note that in calculation 
the effective ionosphere scale height is smoothed over the nearest grid cell and its 
specific values is not crucial. However, for charge exchange between the ions and atoms 
we use the analytical formula for the density of the atmosphere without any smoothing. 
The effective dimensionless diffusion length of the upstream and ionosphere plasmas 
is ld,up = 0.00125, whereas for Io’s body the diffusion length is = 0.0125, where 


Id = 1 /Re and the magnetic Reynolds number is Re = 4nU Q Lo e f f/c 2 . The average 
Lunquist number for the ionosphere may be estimated as 
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Re 


M A 


where the characteristic diffusion and Alfven times are = AnoL 2 /c 2 and ta = L/v a- 
Note that Linker et al. (1998) used the following Lundquist number, Lq = 500 — 1000, 
magnetic Reynolds number, Re = 200 — 400 and diffusion length, l ( j = 0.00125 for 
the background plasma. This value of the diffusion length corresponds to an effective 
conductivity of about 1.2 x 10 8 f2 -1 m -1 far in upstream. Since many plasma and 
atmosphere parameters are still uncertain we have to study a wide spectrum of simple 
models in order to choose the best one for interpretation of the observational data. We 
can present here only a sample of the wide spectrum simulation results which are in 
agreement or in disagreement with observations in order to illustrate the dependence of 
the plasma environment near Io on the input parameters of the problem. The global 
structure of Io’s environment is determined by a set of dimensionless independent 
parameters such as M\, /3 P , /3 e , Mp\/M p , ion production and charge exchange rates, 
diffusion lengths, and the ion gyroradius e = p C i/L. For real values of the magnetic field 
the value of the ion gyroradius is about 8 km which is calculated by the use of the 
local bulk velocity. The dimensionless ion gyroradius and grid spacing have the values 
e = 0.0045 and A X /L = 0.1. In order to study the ion kinetic effect we have to resolve 
the ion gyroradius on the grid. For this reason we use the artificially increased value of 
this parameter, e = 0.126, that is about of a value of A X /L now. Although, such way 
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allows us to study the some kinetic effects for realistic physics we have to extrapolate 
our results into a realistic scale. By scaling the gyroradius in such a way, we preserve 
the ratio of M A , Ms, /3 p , /3 e , and more importantly accurately preserve any anisotropy 
of the ion distribution function with respect to the magnetic field. 

3.1. Global Structure of the Io’s Environment 

Let us consider first the global picture of the interaction of the plasma torus with Io 
in a case with an ionization rate <3j on = 3.03 x 10 27 s -1 (run ceie2, Table 1) and charge 
exchange of form ^ + B exp . The maximum value of the atmosphere density 
is n atm os = 5 x 10 7 cm~ 3 and diffusion lengths are Zd.up = 0.00125 and Z^io = 0.00125. 
Figures 2 and 3 demonstrate the asymmetrical distribution of the torus plasma ion (top) 
and pickup ion density (bottom) in the x-z and y-z planes. One can see the increase 
of the incoming ion density upstream of Io. The pickup ion motion is determined 
mainly by the electromagnetic drift. The motion along the magnetic field is due to 
the thermal velocity and the gradient of the electron pressure. The asymmetrical 
distribution of the incoming ions in the y-z plane may be explained by an existence 
of the B z component of the upstream magnetic field. The inclination of the magnetic 
field results in the asymmetrical boundary condition for ion dynamics (penetration and 
reflection) in the Io’s ionosphere and the asymmetrical Alfven wing. The incoming 
ions flow around the Alfven wing so that only a small portion of the incoming flux 
penetrates the wing, resulting in a decrease in the incoming ion density. This effect is 
stronger in upper half-plane (y > 0) (Fig. 3, top) because the Alfven wing has a stronger 
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front transition in this region than in a lower half-plane. The pickup ion distribution, 
Fig. 3 (bottom), gives the correct value for the inclination of the Alfven wing, 21-21.5 
degrees. The density profile is a little bit disturbed near the side boundaries, however, 
this perturbation does not affect the region close to Io as discussed already. Figures 4 
and 5 show the distribution of the electric and magnetic fields. The asymmetry of the 
distributions in E and B appears to be caused by finite gyroradius effects of incoming 
and pickup ions. A weak perturbation of the magnetic field was observed near the 
ionosphere of Io: compression of the magnetic field upstream and decompression in the 
plasma wake. 

Figure 5 also shows the formation of a strong Alfven (and whistler) wing in the 
direction of the magnetic field. The perturbation of the electric field inside the wings is 
very strong and it may affect ion dynamics so that particles flow around the wings. The 
formation of the Alfven wing in a subalfvenic flow near Io was studied first analytically 
by Neubauer (1980). A excitation of a whistler wave near a plasma cloud was studied by 
using 3-D hybrid simulation in Lipatov (2002). The above wave propagation is closely 
connected with the generation of low frequency waves by the harmonic dipole (local 
source) in the magnetized plasma. The first analytical studies of these effects may be 
found, for example in (Van’yan and Lipatov, 1972, and references therein). 

Figures 6 and 7 show the velocity arrows of incoming and pickup ion velocities. The 
incoming ions flow around the effective obstacle that is produced by pickup ions and 
the ionosphere. The pickup ions flow from the “corona” across the magnetic field due to 
electromagnetic drift whereas the motion along the magnetic field is determined by the 
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thermal velocity of ions and the electron pressure. Figure 8 shows the two-dimensional 
cross section for total ion density in the x-y plane. One can see the asymmetry of 
distribution relative to the x-axis due to the angle between the bulk velocity and the 
magnetic field upstream and to the y-axis due to effects of the finite ion gyroradius. 

3.2. Effect of Electron Temperature on the Plasma Environment 

Our model takes into account the three characteristic electron temperatures that 
do not play any role in past MHD models: a) the temperature of electrons in the plasma 
torus, Te (Up ; b) the temperature of electrons that are created together with pickup ions, 
T e pi; c) the temperature of electrons that are created with ions in the ionosphere, 

T e< iono- Note that in all cases discussed here the pickup ions are generated only in the 
exosphere, i.e W ext = 0. 

Let us consider the case (a). Figure 9 shows the one-dimensional cuts of the total 
ion density, temperature and magnetic field along the x-axis for z = 1.5 x ri 0 and 
y — 0. One can see two peaks in the distribution of the density, each of which has a 
width of about ri 0 and maximum densities of about 4 and 6.5 relative to the upstream 
density. The depletion of the density at x = 2 L may be explained by the following. 
When the supersonic flow passes around the obstacle, a wake with a decreased density 
is formed. The simple example of a such void is presented by the lunar plasma wake. 
However, in our case an asymmetrical pickup ion high density obstacle provides an 
asymmetrical void in Io’s plasma wake. The temperature profile has two peaks with 
maximum temperatures of about 2.3 and 0.7. These peaks do not correlate with the 
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peak in density because the temperature of the ions is determined by the heated ions 
from the incoming flow and pickup ions. The magnetic field profile shows a decrease 
in the magnetic field with a minimum value of about 0.5 in the plasma wake and with 
some significant oscillations (Fig. 9). Figure 10 shows a two-dimensional cross section of 
the total ion density across the plasma wake ( z = 1.5 x ri 0 ). One can see three strong 
maxima. This distribution is determined by finite gyroradius effects of ions. 

If we take into account the temperature of electrons that are created with pickup 
ions (case (b)) the distribution of the ions in the plasma wake may be changed 
significantly. The one-dimensional total density profile has one peak with a value of 
about 8 (Fig. 11). The temperature profile has also one strong peak with a value of 
about 1.9 and the magnetic field has a depletion with a minimum value 0.7 (Fig. 11). 
Figure 12 shows a two-dimensional cross section of the total density across the plasma 
wake. One can see a strong peak narrow in x-direction and wide in y-direction. 

In case (c), when the temperature of electrons inside the ionosphere is also taken 
into account, we have the following distribution of plasma parameters in the wake. 
Figure 13 shows one-dimensional cuts of the total density, temperature and magnetic 
field. The density profile has a wide peak with a thickness of about 4ri 0 and a peak 
value of n rs 6. The temperature profile has a two-peak distribution. The value of 
the temperature at the two peaks is about 0.6 and 0.55. The peaks are separated by 
a distance of 4 x ri 0 . The magnetic field profile has a minimum at the plasma wake 
with B = 0.55 (Fig. 13). Figure 14 shows the two-dimensional cross section of the 
total density across the plasma wake. One can see a strong peak which is narrow in 
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the y-direction and wide in the ^-direction. By introducing cool electrons inside the 
ionosphere we are trying to account in an appropriate way for the cooling effect on 
electrons through electron-collisions close to Io. 

Finally, we show the two-dimensional cuts of pickup ion density for all three cases, 
Fig. 15 (top, middle, bottom). One can see that the temperature of electrons which 
are created in different regions of Io’s environment may affect strongly the pickup ion 
distribution. Figure 15 demonstrates the importance of the electron pressure in the 
region close to the Io. 

In the case of the absence of pressure of pickup and ionospheric electrons, the 
interaction of the plasma torus ions with Io is strongly asymmetrical in the x-z plane 
due to the finite ion gyroradius. In the plasma wake a ’’space-mixing” of the pickup ions 
is observed that results in the formation of an asymmetrical tongue-type distribution 
in the x-z plane. The effusion of pickup ions in the y direction is also weak because of 
small values of A p e and the velocity of the pickup ions in y direction, Fig. 15, (top). 

In the opposite case, when the gradient of the pressure of the pickup electrons and 
the pressure of ionospheric electrons cause a strong electrostatic field, the pickup ion 
density distribution in the x-z plane becomes more symmetrical and wider across the 
wake compared with the above case (cf. Fig. 15 (top) and (bottom)). The gradient of 
the electron pressure also causes a strong effusion of the pickup ions along the Alfven 
wing (see, Fig. 15 (bottom) and (top)). 

Finally, in the intermediate case when the pressure of ionospheric electrons is small 
the asymmetrical distribution of the pickup ion density with a strong effusion along the 
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Alfven wing is observed, Fig. 15 (middle). So, we can conclude that the pressure of the 
pickup ions is responsible for effusion processes while the pressure of the ionospheric 
electrons is responsible for the formation of the symmetrical distribution with the 
maximum in density located near the (y=0) plane. 

All three cases, (a), (b) and (c), demonstrate strong asymmetries in the distributions 
of pickup ion density, temperature and magnetic field in comparison with MHD and 
electrodynamic models. 

3.3. Effect of Pickup Ion Injection Distribution on the Plasma Environment 

To study the role of pickup ions from the halo we simulated Io’s environment 
for cases when the ion production rate in the halo corresponds to W ext — 10% and 
W ext = 30% of the total (halo plus exosphere) ion production rate. Note that all other 
parameters are the same in these cases. 

Figure 16 (solid line) presents one-dimensional cuts for density, temperature and 
magnetic field for W ext = 10% and /3 e ,iono = 0.25/? e . The density profile has a maximum 
value of about n « 7.2 with a characteristic width of the peak of 3 x ri 0 , that is higher 
than the observational value, n as 5.1. The temperature profile shows a two-peak 
distribution. The maximum values in the peaks are T/T 0 rs 3 and T /Tq rs 2.1 that 
correspond to the observational data (Frank et al., 1996; Kivelson et ah, 1996). The 
distance between peaks is about 3 x ri 0 that is a little bit wider than observed value, 

2 x ri 0 . The magnetic field profile has a minimum value in the plasma wake of B « 0.5 
that also corresponds the minimum value of the B y in the observational data. Figure 16, 
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(dotted line) shows the one-dimensional cuts for density, temperature and magnetic 
field for W ext = 30% and fi e ,i ono = 0.25/3 e - The density profile has a maximum value 
of about n « 5.83 with a characteristic width of the peak of 4 x ri 0 , which is a little 
bit wider than that in Fig. 16, (solid line) and is a little higher than observed. The 
temperature profile also shows a two-peak distribution. The maximum values in the 
peaks are T/To ~ 6.8 and T /T 0 ~ 2.4. Note that the first maximum in the temperature 
is higher than observed (Frank et al., 1996; Kivelson et ah, 1996). The magnetic field 
profile shows a minimum value of the magnetic field in the plasma wake B « 0.56 that 
corresponds the observational data. 

3.4. Effect of Ionospheric Conductivity on the Plasma Environment 

In the previous cases we have assumed a high conductivity for the ionosphere and 
Io’s body. The realistic models of Io may include a conducting core, surrounded by a 
poorly conducting mantle (a Moon-like model), or Io may be considered as a poorly 
conducting body. In this section we shall model the Io as a poorly conducting body. 
Note that all other parameters in this case are the same as that shown in Fig. 16, (solid 
line) except the higher diffusion length. Figure 16 (dashed line) shows one-dimensional 
cuts for the density, temperature and magnetic field for the case with a diffusion length, 
ld,i<mo — 0.025. Note that the lines in Fig. 16 corresponds the simulation for 
different sizes of the computational domain. The value of the peak density is 
about n « 7.5 with the same thickness of the peak, « 3r Io as in Fig. 16, (solid line). 
The maximum values of the peak temperature, 3 and 2 are approximately the same as 
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in Fig. 16, (solid line) (see Fig. 16, (dashed line)). However, the magnetic field profile 
shows much stronger field variation in the external region of the plasma wake than in 
Fig. 16, (solid line). The analysis of the magnetic field inside Io shows that the reduced 
conductivity (Zd,iono = 0.025, Fig.18 (dashed line)) near ionosphere and inside Io may 
reduce the magnetic field inside Io by 20 — 40% and the peak plasma density at the 
Galileo trajectory by 5% in comparison with the model with Zd.iono = 0.0025). 

3.5. Comparison with Galileo Observational Data 

The result of the measurements by the particle and field instruments on the Galileo 
Orbiter during the December 1995 flyby of Io provided new and important information 
with which realistic simulations for the plasma interaction can be tested. Along that 
trajectory physical signatures of the wake were seen as a broad depression in the 
magnetic field (Kivelson et al., 1996), sharp peaks in the ion (Frank et al., 1996) and 
electron (Gurnett et al., 1996) densities, a slowing of the plasma in the core of the wake, 
a deep ion temperature decrease in the center of the wake, and a large (factor of 3) 
ion temperature rise in the flanks of the wake (Frank et al., 1996). The magnetic field 
perturbation was broader spatially than the density peak, and showed a double-reversed 
structure, whereby the perturbation (as defined by the difference from the outer Jovian 
B-field value) was actually weaker right near the close-approach point than it was 
somewhat adjacent to the center of the wake. For comparison of the our computational 
model with observation data we made a run with the following plasma parameters: 
/3e,pi == 0.125/3 e>U p, /3e,iono — /5e,up> u.j 0no = 100rio, -Hatmos = 0.06ri o , W ex t — 5%, and 
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Z d>up = 0.00125. The conductivity of Io is chosen so that Zd,io = 0.0125. 

Let us consider first a case in the absence of charge exchange processes (case nee 
Table 1). The total ion production rate was Qion = 4.8 x 10 27 s _1 . Fig. 17 shows the 
comparison of plasma parameters and magnetic field components from the simulation 
model and observation along the Galileo trajectory (Pass 10). Note that we used the 
interpolation of the grid values of the plasma parameters and the magnetic 
field into the position of Galileo spacecraft. We can see that the density peak 
is slightly shifted to the right in comparison with observed one. The left peak in the 
temperature profile is a little bit narrower than the observed one while the right peak 
is a little bit higher than the observed one. The simulation yields a little bit smaller 
value of By at large x, x >> r Io , due to perturbation of the electromagnetic field in the 
region above ( y > 0) and below ( y < 0) the equatorial plane. Unlike the MHD models 
the hybrid model yields a B y profile with reverse structure in the middle of the plasma 
wake. This type of behavior may be explained by the diamagnetic and accelerational 
drift currents in the plasma wake, which are modeled naturally in our kinetic description 
for ions. The maximum value of these currents is located near the equatorial plane 
in the boundary layer (interface) between the external plasma flow and pickup ions 
in the plasma wake. The B z profiles in these simulations correspond very well to the 
observed data however the left maximum in the magnetic field is smaller than observed 
one. Unfortunately fluctuations in magnetic field components are not very small due to 
“shot” noise. So the lack in agreement between observed and computed values of B x is 
probably just the result from less than perfect plasma simulation model (“shot” noise), 
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and the simplified models of the ionosphere and Io’s body. 

Let us consider now the cases with a charge exchange rate that corresponds to 
the maximum value of the density of the atmosphere, n atmos = 10 8 cm -3 . We assume 
here that the charge exchange process is due to only the exponential part of the neutral 
atmosphere. 

Fig. 18 (solid line) shows the comparison of plasma parameters and magnetic field 
components from our simulation model and the observations in the case with a total ion 
production rate, Q j on = 2.27 x 10 27 s -1 . In this case the density profile has a smaller 
effective value but the total temperature is a little bit higher. The left peak in the 
temperature profile is also a little bit narrower than in the observational data while the 
right peak is a little bit higher than the observed one. The perturbation in B y is a little 
bit smaller than in the case of absence of charge exchange, (cf. Fig. 18 (solid line) and 
Fig. 17). 

In case of higher total ion production rate, Q\ OD = 2.57 x 10 27 s -1 (Fig. 18, (dotted 
line)) the density profile has a higher effective value but the total temperature is a little 
bit lower. We can see that the density profile has a depletion on the left side and the 
width of the peak is smaller than in the case with no charge exchange, (cf. Fig. 18, 
(dotted line) and Fig. 17). The left peak in the temperature profile is also a little bit 
narrower than in the observational data while the right peak is a little bit higher than 
the observed one. The perturbation in B y is a little bit smaller than in case of absence 
of charge exchange (cf. Fig. 18, (dotted line) and Fig. 17). 

In a case with a higher charge exchange rate that corresponds to the maximum 
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value of the density of atmosphere, n atmos = 5 x 10 8 cm -3 and a total ion production 
rate, Q[ 0n = 2.1 x 10 27 s _1 we have a satisfactory agreement in the total density profile, 
a smaller total temperature and no agreement in the magnetic field component B y . 

Let us consider the cases with a high charge exchange rate that corresponds to the 
maximum value of the density of atmosphere, n atmos = 10 9 cm~ 3 and an ionization rate 
of Qion = (1-4 - 2.1) x 10 27 s~L In case of low ionization rate, Q ion = 1.4 x 10 27 s -1 , 
the density peak is higher and thinner than the observed one whereas the temperature 
profile is in a good agreement with the observation. The magnetic field B y profile is 
only in qualitative agreement with observation. In the case of a higher ionization rate, 
Qion — 2.1 x 10 27 s~ 1 , the density has a strong double peak structure but the temperature 
profile is in agreement with the observations. However, the magnetic field profile has no 
agreement with observational data. 

In previous cases we assumed that the charge exchange process is only due 
to the lower altitude, exponential part of the neutral atmosphere. Let us consider 
now the results of simulation that includes also the charge exchange due to the 
extended atmosphere which is distributed as r~ 2 . Figure 19 shows the profiles of 
the total density, temperature and magnetic field for different maximum values of 
the density of atmosphere, n atm os = (5 x 10 7 — 10 9 )cm" 3 , and an ionization rate, 

Qion = (3.03- 4.04) x 10 27 s~L 

Let us consider first the cases with a low maximum value of the density of 
atmosphere, n atm os = 5xl0 7 cm“ 3 , and an ionization rate ofQ ion = (3.03-4.04) xlO 27 s' 1 . 
In case of lower ionization rate, Q ion = 3.03 x 10 27 s~ 3 , the density profile is in agreement 
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with the observation, but the temperature is a little bit higher in the peaks (Fig. 19, solid 
line). The magnetic field profile ( B y ) is in a qualitative agreement with the observation 
(Fig. 19, solid line). In the case of higher ionization rate, Q ion = 4.04 x 10 27 s _1 , the 
density profile is a little bit higher in peak than the observed one and the temperature 
is a little bit lower in the left peak. The magnetic field profile ( B y ) is much smoother 
than in the observation (not shown in Fig. 19). 

In the cases with a higher maximum atmosphere density, n atmos = 10 8 cm -3 , the 
maximum value of the wake density becomes much higher than in the observation. 
Although the temperature profiles are in good agreement with the observations, the 
magnetic field component B y is much smoother in the case of Q ion = 3.03 x 10 27 s _1 
(Fig. 19, (dotted line)), and it is in qualitative agreement with the observations in the 
case with Q; on = 4.04 x 10 27 s _1 (not shown in Fig. 19). The increase in the atmosphere 
maximum to n atmos = 5x10® s -3 results in an increase of the density to more than 
twice the observed value. In this case we have Qj on = 4.04 x 10 27 s _1 with a total charge 
exchange rate of Q exch = 5.8 x 10 28 s _1 , Fig. 19, (dashed line). 

In case of a high maximum atmosphere density, n atmos = 10 9 cm -3 there is 
no similarity between simulation and observation. We simulated the cases with 
ionization rates, Qion = (2.1 — 4.04) x 10 27 s _1 . In the case of low ionization rate, 

Qion = 2.1 x 10 27 s -1 , the density profile has two thin peaks (Fig. 19, (dot-dashed 
line)) while the magnetic field profile ( B y ) has a strong variation (Fig. 19, (dot-dashed 
line)) that does not correspond the observation. In case of moderate ionization rate, 
Qion = 3.03 x 10 27 s _1 , the density profile has one peak with a value, ntotai > 10 uq 
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that is much higher than observed one whereas the magnetic field profile ( B y ) is much 
smoother than observed one. In the case of a high ionization rate, Qion = 4.2 x 10 27 s _1 , 
the density profile has one peak with a value, n to tai ~ (8 — 10)n 0 that is much higher 
than observed whereas the magnetic field profile (B y ) is much smoother than observed. 

One of the main issues of the plasma torus - Io interaction is the question about 
the real obstacle in this interaction. There are two possible candidates for this problem 
- (1) mass loading by the pickup ions originally produced from Io’s atmosphere by the 
ionization processes or charge exchange processes between the incoming plasma torus 
ions and the neutral atmosphere, or (2) charge exchange processes between the pickup 
ions originally produced from Io’s atmosphere and the neutral atmosphere. Table 
1 gives a summary of results of simulations for different models that are 
considered in this section. 

In the case of the models with a high density of Io’s atmosphere (10 9 cm -3 ) the 
charge exchange rate becomes very high, Qexch ~ (6.7 — 9) x 10 28 s _1 (cases ceie6, 
ceie7 and seie8, Table 1). Note that in these cases the charge exchange rate for 
incoming ions is much smaller than for pickup ions. We compared the charge exchange 
rate in the external (r > 4ff atmos ) and the internal (r < 4 x 77 atmos ) region. The charge 
exchange rate for incoming ions in the external region has approximately the same value 
as a charge exchange rate in the internal region. However, the charge exchange rate for 
pickup ions in the external region is approximately 20 times smaller than the charge 
exchange rate in the internal region. 

The magnetic field profile ( B y ) does not match the observational data very well 
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for a small ionization rate, Qj on = 2.1 x 10 27 s -1 (case ceie3, Table 1), whereas the 
total density value in the plasma wake is much more than observed one in the case of 
higher ionization rate, Q i 0n = (3.03 - 4.04) x 10 27 s~ x (cases ceiel and ceie2, Table 
1). Hence, the models with a high charge exchange rate, Q exC h = (6.5 - 9) x 10 28 s _1 , 
are not realistic because they cannot explain the observational data. 

The model with n atmos = 5 x 10 8 cm -3 (case ceie5, Table 1), gives the total charge 
exchange rate Qexch = 5.8 x 10 28 s _1 . The models with n atm0 s = 10 8 cm -3 (cases ceie3 
and ceie4, Table 1), give total charge exchange rates of Q ex ch = (1 09- 1.35) x 10 28 s _1 
And, finally, the models with n atmos = 5 x 10 7 cm~ 3 , give total charge exchange rates of 
Qexch = (6.12 - 6.54) x 10 27 s _1 (cases ceiel and ceie2, Table 1). 

The analysis of the models shows that two cases that have a good fit for 
observational data. In the model without charge exchange (case nee, Table 1) the 
simulation results describe well enough the main observational data - a total density, a 
total temperature and the magnetic field ( B y ) (Fig. 17). The more realistic model with 
n atmos = 5 x 10 7 cm -3 , and ionization rate Qj on = 3.03 x 10 27 s _1 (case ceiel, Table 1), 
also fits satisfactorily the observational data (Fig. 19, solid line). This model gives the 
charge exchange rates Q p , e xch = 3.16 x 10 27 s -1 and Q c ,exch = 3.13 x 10 27 s _1 for plasma 
torus ions and pickup ions respectively and would correspond to an MHD model with 
total fresh ion mass-loading rate and the charge exchange rate, which contributes to the 
momentum and energy friction terms, of 1.2 x 10 28 s _1 (Combi, Gombosi and Kabin, 
2002). So, we find that the interaction is dominated in roughly equal parts between 
primary torus ion charge exchange and secondary pickup ion charge exchange. 
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4. Conclusions 

3D Boltzmann simulations of the interaction of the Jovian plasma torus with Io, 
have demonstrated several new features: 

• The whistler (lead front) and quasi-stationary Alfven waves (wings) strongly affect 
plasma flow around Io. 

• The effect of the finite ion gyroradius that results in formation of the asymmetrical 
boundary layer in the vicinity of Io’s ionosphere is important. The plasma 
parameters have a strongly asymmetrical distribution across Io’s wake. The 
kinetic behavior of ion dynamics reproduces the inverse structure of the magnetic 
field (due to drift current) which cannot be explained by standard MHD or 
electrodynamic simulations which do not account for anisotropic ion pressure. 
The diamagnetic effect of non-isotropic gyrating pickup ions broadens the 
B-field perturbation and produces increased temperatures in the flanks of the 
wake, as observed by the Galileo spacecraft, but not explained by previous 
simulations. Note that two-fluid simulation (Saur et al., 2002) produces 
the double-peak signature with spatial scale much smaller than in 
observation. 

• The cold, dense wake is produced, as in MHD, but was not produced in 
electromagnetic simulations without ad hoc addition of bidirectional electrons. 


• The values of temperatures of the electrons which are created and cooled by 
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collision with neutrals in the exosphere and inside the ionosphere may strongly 
affect the pickup ion dynamics along the magnetic field and consequently the pick 
up distribution across the wake. In the absence of an observed global picture of the 
plasma distribution, the simulation serves to demonstrate the wide range of global 
configurations that are possible for various electron temperature descriptions. A 
full hybrid simulation for the plasma including an electron temperature treatment 
that accounts for electron-neutral collision and a detailed neutral description 
would be required to produce an accurate global picture from the first principles. 

• The effective conductivity of Io’s ionosphere may change strongly the distribution 
of the magnetic field near and inside Io. The reduced conductivity (i l dif = 0.025) 
near the ionosphere and inside the Io may reduce the magnetic field inside Io by 
20 — 40% and the peak of the plasma density at the Galileo trajectory by 15%. 

• The best models that fit the observational data well are: the model without charge 
exchange and ionization rate Qi 0n = 4.8 x 10 27 s -1 , and the model with ionization 
rate Q ion = 3.03 x 10 27 s _1 and the total charge exchange rate Q e xch = 6.3 x 10 27 s -1 . 
The best MHD model had a total fresh ion mass-loading rate of 6 x 10 27 s -1 and 
a charge exchange rate, which contributes to the momentum and energy friction 
terms, of 1.2 x 10 28 s -1 (Combi, Gombosi and Kabin, 2002). The reason for the 
factor of 2 difference results from the above-described pressure anisotropy effect 
in the kinetic simulation, and the resulting tighter plasma distribution near the 


equator plane defined by Io. 
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• In the presented calculations the value of the ion gyroradius was much 
larger than the realistic value. So we have to investigate the effects of 
realistic values of the ion gyroradius in future simulations. We expect 
that the smaller ion gyroradius may result in smaller asymmetry of 
the global picture. However, we do not expect the strong changes in 
the plasma and the magnetic field profiles because the used value of 
the ion gyroradius is much smaller than the characteristic scale of the 
problem - the radius of Io. These profiles are controlled primary by the 
inter-penetration of the torus plasma and the pickup ions. 
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5. Figure captions 

Figure 1 . Galileo trajectory close to Io and the system of coordinates. 

Figure 2. Incoming (top) and pickup (bottom) ion density in the x-z plane. The 
case with Q ion = 3.03 x 10 27 s _1 , n atmos = 5 x 10 7 cm~ 3 , /3 e<P[ = /3 e , /3 e>iono = 0.25 /3 e , 
W ext = 5%, and H atmos — 0.06. See explanation in Fig. 6. 

Figure 3. Incoming (top) and pickup (bottom) ion density in the y-z plane for the 
same parameters as Fig. 2. See explanation in Fig. 7. 

Figure 4. Electric (top) and magnetic (bottom) fields in the x-z plane for the 
same parameters as Fig. 2. Figure shows a strong asymmetry of the electromagnetic 
field due to finite gyroradius effects. 

Figure 5. Electric (top) and magnetic (bottom) fields in the y-z plane for the same 
parameters as Fig. 2. Figure shows the formation of an Alfven wing in the direction of 
the main magnetic field. 

Figure 6. Incoming (top) and pickup (bottom) ion velocity arrows in the x-z plane 
for the same parameters as Fig. 2. The figure demonstrates a flow of pickup ions from 
the “corona” across the magnetic field. The incoming ions flow around the effective 
obstacle that is produced by pickup ions and ionosphere. 

Figure 7. Incoming (top) and pickup (bottom) ion velocity arrows in the y-z 
plane for the same parameters as Fig. 2. The figure demonstrates a strong expansion of 
pickup ion “corona” along the magnetic field line. The incoming ions flow around the 


region of extended “corona” . 
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Figure 8. Two-dimensional section for total density in the x-y plane for the same 
parameters as Fig. 2. 

Figure 9. Total density, temperature and magnetic field along the x-axis. The 
case with Pe,pi — 0, Pe,i ono ~~ IF ext ~ 0%, and H a tmos — 0.06. 

Figure 10. Two-dimensional section for total density in the x-y plane. The case 
with Pe iPI = 0, Pe^ono = 0, W ext = 0%, and H atmos = 0.06. 

Figure 11. Total density, temperature and magnetic field along the x-axis. The 
case with Pe,pi — Pet Pe,iono — 0 , W ex t — 0%, and H a tmos — 0.06. 

Figure 12. Two-dimensional section for total density in the x-y plane. The case 
with Pe,pi = Pet Pe,iono = 0, W ext = 0%, and H atmos = 0.06. 

Figure 13. Total density, temperature and magnetic field along the x-axis. The 
Case with Pe,PI = Pet Pe,iono = Pet W ex t = 0%, H a tmos ~ 0.06. 

Figure 14. Two-dimensional section for total density in the x-y plane. The case 
with Pe,PI — Pet Pe,iono = Pet ^ext = 0%, Hatmos ~ 0.06. 

Figure 15. Two-dimensional sections for pickup ion density in the x-y plane. 
The cases with (top) p giPI = 0, p e ,imo = 0; (middle) p e<PI = P e , p e , iono — 0; (bottom) 
Pe,pi = Pet Pe.iono = 0.25 P e . W ext = 0%, and H atmos = 0.06. Figure demonstrates 
asymmetry of pickup ion density across the wake. 

Figure 16. Total density, temperature and magnetic field along the x-axis for 
Pe,pi = Pet Pe,i<mo = 0.25/5 e , H atmos = 0.06, l d ,u P = 0.0025. (solid line) W ext = 10%, 
ld,iono = 0.0025; (dotted line) W ext = 30%, ld,iono = 0.0025; (dashed line) W ext — 10%, 
ld,iono — 0.025. 
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Figure 17. Comparison of Galileo PLS (10 pass) [ Frank et al, 1996] and MAG 
[ Kivelson et al., 1996] (the open circles) data with the Io hybrid model results in absence 
of charge exchange. Q = 4.8 x 10 27 s -1 , /3 e ,p/ = 0.125 /3 e , f3 e ,iono = 0.125/? e , W ext = 5%, 
Hatmos = 0 06, and n iono = 800(100)n o . Id, up = 0.00125 and ldj 0 = 0.0125. 

Figure 18. Same as Fig. 17, but with charge exchange (Sex P^~^j), 
n atmos — 10 8 cm -3 : (solid line) Q = 2.27 x 10 27 s _1 ; (dashed line) Q — 2.57 x 10 27 s _1 . 

Figure 19. Same as Fig. 17, but with charge exchange (^ + Sexp ) : (solid 
line) n atmos = 5 x 10 7 cm~ 3 , n iono = n 0 , Q = 3.03 x 10 27 s _1 , Q v , exc h = 3.169 x 10 27 s -1 , 
Qc,exch = 3.13 x 10 27 s -1 , Idjo = 0.00125; (dotted line) n atmos = 10 8 cm“ 3 , 
niono = n 0 , Q = 3.03 x 10 27 s _1 , Q v , exc h = 5-27 x 10 27 s _1 , Q c ,exch ~ 5.6 x 10 27 s -1 , 
h,io = 0.00125; (dashed line) n atmos = 5 x 10 8 cm -3 , ni ono — n 0 , Q = 4.04 x 10 27 s _1 , 
Qp,exch = 1-83 x 10 28 s _1 , Q c ,exch = 3.97 x 10 28 s _1 , Idjo = 0.0125; (dot-dashed 
line) n a t mos = 10 cm , n^ono — no, Q ~ 2.1 x 10 s , Q p ,exch — 2.95 x 10 s , 
Qc,exch = 3.75 x 10 28 s _1 , l d jo = 0.0125. 



Table 1 . Dependence of torus ion charge exchange rate and pickup ion charge exchange rate 
versus the maximum value of atmospherical density and ionization rate (Q 0 = 10 27 s -1 ). 


* denote the regimes that fit satisfactory the observational data 

CaSe lQ^cm- 3 Qp 4o Ch Qc $o Ch Comments 

model without charge exchange 

nee 4.8 0 0 a density is a little bit wide and it is in good agreement with an observation, 

temperature and magnetic field are in a good agreement with an observation, Fig. 17* 
models with charge exchange in the internal region (neutral ~ B exp Hfl ) 

“atmos 

ceil 1 2.27 density, temperature and B y in satisfactory agreement with observation, Fig. 18, 

(solid line) 

cei2 1 2.57 density and temperature in satisfactory agreement with observation, B y component 

is a weaker than observed, Fig. 18, (dashed line) 

cei3 5 2.1 density in good agreement with observation; temperature and B y are not in agreement 

with observation 

cei4 10 1.4 12.6 14 density with strong narrow peak, temperature in satisfactory agreement with 

observation; B y is weaker than observed 

cei5 10 2.1 12 20.2 density with two peaks, temperature is in satisfactory agreement with observation; 

By has much stronger depletion than in observation 

models with charge exchange in the internal and external region (n neutra i k, 4 + Bex p r u ~ r k » ) 

r ** atm os 

ceiel 0.5 3.03 3.169 3.13 density, temperature and B y in good agreement with observation. Fig. 19*, (solid line) 

ceie2 0.5 4.04 2.67 3.87 density and temperature in satisfactory agreement with observation, B y is a weaker 

than observed 

ceie3 1 3.03 5.27 5.6 density is higher than the observation, variation in B y is a weaker than observed, 

Fig. 19, (dotted line) 

ceie4 1 4.04 5.55 7.93 density is higher than observation, variation in B y is weaker than observed 

ceie5 5 4.04 18.3 39.7 density is much higher than in an observation, variation in the B y is strong but profile 

is wider than observed, Fig. 19, (dashed line) 

ceie6 10 2.1 29.5 37.5 two-peak density profile, variation in the B y is much stronger than observed, 

Fig. 19, (dot-dashed line) 

ceie7 10 3.03 28. 44.8 density is much higher than in observation, variation in the B y is much weaker than 

observed 

ceie8 10 4.04 27.5 60. density is much higher than in observation, variation in B y is much weaker than 


observed 
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Figure 1 . Galileo trajectory close to Io and the system of coordinates. 





Figure 2. Incoming (top) and pickup (bottom) ion density in the x-z plane. The case 
ceiel (Table 1) with /9 e .pi = A>,iono = 0.25 /3 e , W ext = 5%, and I/ atmos = 0.06. See 
explanation in Fig. 6. 
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Figure 3. Incoming (top) and pickup (bottom) ion density in the y-z plane for the same 
parameters as in Fig. 2. See explanation in Fig. 7. 
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Figure 4. Electric (top) and magnetic (bottom) fields in the x-z plane for the same 
parameters as in Fig. 2. Figure shows a strong asymmetry of the electromagnetic field 
due to finite gyroradius effects. 








Figure 6. Incoming (top) and pickup (bottom) ion velocity arrows in the x-z plane for 
the same parameters as in Fig. 2. The figure demonstrates a flow of pickup ions from the 
“corona” across the magnetic field. The incoming ions flow around the effective obstacle 
that is produced by pickup ions and ionosphere. 





Figure 7. Incoming (top) and pickup (bottom) ion velocity arrows in the y-z plane for 
the same parameters as in Fig. 2. The figure demonstrates a strong expansion of pickup 
ion “corona” along the magnetic field line. The incoming ions flow around the region of 


extended “corona”. 
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Figure 8. Two-dimensional cross section for total density in the x-y plane at z = 1.5rj 0 
for the same parameters as in Fig. 2. 
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Figure 9. Total density, temperature and magnetic field along the x-axis at z = 1.5rj 0 . 
The case with /? ejP i = 0, /3 e ,iono = 0, W ext = 0%, and i^atmos = 0.06. 
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Figure 10. Two-dimensional cross section for total density in the x-y plane at z = 1.5ri 0 . 
The case with /3 e>P i = 0, /3 e , i ono = 0, W ext = 0%, and // atm0 s = 0.06. 
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Figure 11. Total density, temperature and magnetic field along the x-axis at z = 1.5ri 0 . 
The case with &,pi = ft, ft,k>no = 0, W^t = 0%, and tf atmos = 0.06. 


N 



X/L 

Figure 12. Two-dimensional cross section for total density in the x-y plane at z = 1.5ri 0 . 
The case with /J e ,Pi = ft, ft,iono = 0, W ext = 0%, and // a tmos = 0.06. 
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Figure 13. Total density, temperature and magnetic field along the rr-axis at z — 1.5ri 0 . 
The case with y3 e ,pi = Pe, Pe,iono = Pe, W ext = 0%, # a tmos = 0.06. 
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Figure 14. Two-dimensional cross section for total density in the x-y plane at z = 1.5ri 0 . 
The case with /3 e ,pi = Pe, Pe,\ono = Pe, W ex t = 0%, H itmos = 0.06. 






Figure 15. Two-dimensional cross sections for pickup ion density in the x-y plane at 
2 = 1.5r Io . The cases with (top) &,pi = 0, /? e ,iono = 0; (middle) &.pi = /3 e ,iono = 0; 
(bottom) /? e ,pi = p e , Pejono = 0.25/? e . W ext = 0%, and ifatmos = 0.06. Figure demon- 


strates asymmetry of pickup ion density across the wake. 
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Figure 16. Total density, temperature and magnetic field along the x- axis at, z = 1.5r Io 
for cases with /? e ,Pi = Pe, Pe,iono = 0.25&, H & tm os = 0.06, / djUp = 0.0025. (solid line) 
W ext = 10%, Zd.iono = 0.0025; (dotted line) W ext = 30%, /d.iono = 0.0025; (dashed line) 
W ext = 10%, l d} iono = 0.025. 
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Figure 17. Comparison of Galileo PLS (10 pass) [Frank et al., 1996] and MAG [Kivel- 
son et al., 1996] (the open circles) data with the Io hybrid model results along the 
trajectory in the absence of charge exchange (case nee, Table 1). /3 e ,pi = 0.125/3 e , 
fiejono = 0.125 fi e , W e x t = 5%, # a tmos = 0.06, and niono = 800(100)n 0 . id, up = 0.00125 and 




^d,io = 0.0125. 
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. 17, but with charge exchange ( B exp hJ 1 I s ), n atmos = 10 8 cm 2 3 * : 

>le 1); (dashed line, case cei2, Table 1). 
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Figure 19. Same as Fig. 17, but with charge exchange {£+B exp ): (solid line, case 
ceiel, Table 1) nj 0n0 = no, id,io = 0.00125; (dotted line, case ceie3, Table 1) /d,i 0 = 0.00125; 
(dashed line, case ceie5, Table 1) Zd,io = 0.0125; (dot-dashed line, case ceie6, Table 1) 


^d,io = 0.0125. 



